clear, clc

%%

fields = 1:8;%1:8;
cols = 'L';%'C':'N';
rows = 15;%3:22;

%%

rpath = '/path/to/root/folder';

nModelPath = [rpath filesep 'segModels' filesep 'modelNuclei2.mat'];
cModelPath = [rpath filesep 'segModels' filesep 'modelCells2.mat'];
lModelPath = [rpath filesep 'segModels' filesep 'modelLDs2.mat'];

platePath = 'path/to/plate/images';

cpath = 'path/to/control(background_signal)/images/folder';

%%

load(nModelPath); nmodel = model;
load(cModelPath); cmodel = model;
load(lModelPath); lmodel = model;

%%

load(cpath)
meanBgSQ = meanAllNeg; % average background signal (to be subtracted before quantification)

resultsPath = [platePath '_Results'];
if ~exist(resultsPath,'dir')
    mkdir(resultsPath)
end

rc = cell(1,length(rows)*length(cols));
count = 0;
for iRow = 1:length(rows)
    r = rows(iRow);
    for iCol = 1:length(cols)
        c = cols(iCol);
        count = count+1;
        rc{count} = {r,c};
    end
end

imSize0 = [2048 2048]; % original image size
imSize = [512 512]; % image size used for ML

avgs = cell(1,length(rc));
parfor iRC = 1:length(rc)
% for iRC = 1:length(rc)
    r = rc{iRC}{1};
    c = rc{iRC}{2};
    fprintf('row %d, col %s\n', r, c);
        
    allFieldsQuants = [];
    outImg = [];
    for iField = 1:length(fields)
%         tic
        f = fields(iField);
        pNC = sprintf('%s - %02d(fld %d wv Red - Cy5).tif',c,r,f); % nuclei and cells
        pLD = sprintf('%s - %02d(fld %d wv UV - DAPI).tif',c,r,f); % lipid droplets
        pSQ = sprintf('%s - %02d(fld %d wv Blue - FITC).tif',c,r,f); % signal to quantify

        NC = imresize(imreadGrayscaleDouble([platePath filesep pNC]),imSize);
        LD = imresize(imreadGrayscaleDouble([platePath filesep pLD]),imSize);
        SQ = imreadGrayscaleDouble([platePath filesep pSQ]);
        
        % -----
        % correcting for auto-fluorescence (no bleed-through detected -- see scratchControls.m)
        SQ = SQ-meanBgSQ;
        % -----

        % segment
        [IndLists, SegRGB] = segmentCellsAndLDs(NC,LD,nmodel,cmodel,lmodel,imSize,imSize0);

        % quantify
        if ~isempty(IndLists)
            [quants,QuantRGB] = quantifySignalInCompartments(SQ,IndLists);
            QuantRGB = imresize(QuantRGB,imSize);
            allFieldsQuants = [allFieldsQuants; [f*ones(size(quants,1),1) quants]];
        else
            QuantRGB = zeros(imSize(1),imSize(2),3);
        end
        outImg = [outImg; [SegRGB QuantRGB]];
%         toc
    end
    
    if ~isempty(allFieldsQuants)
        avgAllFieldsQuants = nanmean(allFieldsQuants(:,3:13));
        T = array2table(allFieldsQuants,'VariableNames',{'field','cell_id','int_int_LD','int_int_notLD','area_LD','area_notLD','avg_int_LD','avg_int_notLD','ratio_AvgIntLD_AvgIntNotLD','medLD','med_notLD','iqrLD','iqr_notLD'});
    else
        avgAllFieldsQuants = [];
        T = array2table(nan(1,13),'VariableNames',{'field','cell_id','int_int_LD','int_int_notLD','area_LD','area_notLD','avg_int_LD','avg_int_notLD','ratio_AvgIntLD_AvgIntNotLD','medLD','med_notLD','iqrLD','iqr_notLD'});
%         T = array2table(allFieldsQuants); % empty table
    end
    avgs{iRC} = {r,c,avgAllFieldsQuants};
    
    writetable(T,[resultsPath filesep sprintf('%s - %02d.csv',c,r)]);
    imwrite(outImg,[resultsPath filesep sprintf('%s - %02d.png',c,r)]);
end

avgs2 = cell(length(rc),13);
for iRC = 1:length(rc)
    avgs2{iRC,1} = avgs{iRC}{1};
    avgs2{iRC,2} = avgs{iRC}{2};
    aafq = avgs{iRC}{3};
    if ~isempty(aafq)
        for iQ = 1:length(aafq)
            avgs2{iRC,iQ+2} = aafq(iQ);
        end
    end
end

T = cell2table(avgs2,'VariableNames',{'row','col','int_int_LD','int_int_notLD','area_LD','area_notLD','avg_int_LD','avg_int_notLD','ratio_AvgIntLD_AvgIntNotLD','medLD','med_notLD','iqrLD','iqr_notLD'});
writetable(T,[platePath '_Averages.csv']);
